perm filename SNEW2.F4[DRW,LCS] blob
sn#079381 filedate 1974-12-13 generic text, type T, neo UTF8
00100 SUBROUTINE SS
00200 COMMON X(100),Y(100),N,X1(512),Y1(512),S(100),K
00300 DIMENSION A(2,200),B(2,200),C(2,200),I(200),I1(200),I2(200)
00400 C FIND MAX AND MIN POINTS AND SET SLOPES
00500 ICOUNT=0
00600 710 DO 100 K=2,N-1
00700 IF(Y(K).GT.Y(K-1).AND.Y(K).GT.Y(K+1)) GO TO 20
00800 IF(Y(K).LT.Y(K-1).AND.Y(K).LT.Y(K+1)) GO TO 20
00900 IF(X(K).GT.X(K-1).AND.X(K).GT.X(K+1)) GO TO 30
01000 IF(X(K).LT.X(K-1).AND.X(K).LT.X(K+1)) GO TO 30
01100 GO TO 100
01200 20 S(K)=0
01300 ICOUNT=ICOUNT+1
01400 I(ICOUNT)=K
01500 GO TO 100
01600 30 S(K)=1001
01700 ICOUNT=ICOUNT+1
01800 I(ICOUNT)=K
01900 100 CONTINUE
02000 ICOUNT=ICOUNT+1
02100 I(ICOUNT)=N
02200 C FIND POINTS SUCH THAT X(K)=X(K+1) OR Y(K)=Y(K+1)
02300 DO 900 K=2,N-1
02400 IF(Y(K).EQ.Y(K+1)) I2(K)=1
02500 IF(X(K).EQ.X(K+1)) I2(K)=1
02600 900 CONTINUE
02700 C SET THE SLOPES BETWEEN
02800 IA=1
02900 IX=0
03000 110 IF(IX.EQ.ICOUNT) GO TO 200
03100 IX=IX+1
03200 IA=IA+1
03300 120 IF(IA.EQ.I(IX)) GO TO 110
03400 S(IA)=(Y(IA+1)-Y(IA-1))/(X(IA+1)-X(IA-1))
03500 IA=IA+1
03600 GO TO 120
03700 200 CONTINUE
03800 C SET BEGIN AND END SLOPES
03900 SZ=1002001.
04000 K=1
04100 IF(S(K)**2.EQ.SZ.OR.S(K+1)**2.EQ.SZ) CALL S101
04200 CALL S102(X(1),Y(1),X(2),Y(2),S(1),S(2))
04210 IF(S(K).EQ.-1001.) S(K)=1001.
04220 IF(S(K+1).EQ.-1001.) S(K+1)=1001.
04300 K=N-1
04400 IF(S(K)**2.EQ.SZ.OR.S(K+1)**2.EQ.SZ) CALL S101
04500 CALL S102(X(N),Y(N),X(N-1),Y(N-1),S(N),S(N-1))
04510 IF(S(K).EQ.-1001.) S(K)=1001.
04520 IF(S(K+1).EQ.-1001.) S(K+1)=1001.
04600 C RESET THE SLOPES FOR STRAIGHT LINES
04700 DO 610 K=1,N-2
04800 U1=(Y(K+1)-Y(K))/(X(K+1)-X(K))
04900 U2=(Y(K+2)-Y(K+1))/(X(K+2)-X(K+1))
05000 IF(ABS(U1-U2).GT..0001) GO TO 610
05100 I1(K)=1
05200 I1(K+1)=1
05300 S(K)=U1
05400 S(K+1)=U1
05500 S(K+2)=U1
05600 610 CONTINUE
05700 C ADD POINTS
05800 K1=N-1
05900 DO 300 K=1,K1
06000 IF(I1(K).EQ.1) GO TO 300
06100 IF(I2(K).EQ.1) GO TO 300
06110 FLG=0.
06120 IF(S(K).EQ.0..AND.S(K+1).EQ.0.) FLG=1.
06130 IF(S(K).EQ.1001..AND.S(K+1).EQ.1001.) FLG=2.
06140 IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 205
06200 SZ=1002001.
06300 IF(S(K)**2.EQ.SZ.OR.S(K+1)**2.EQ.SZ) CALL S101
06400 U=(Y(K)-Y(K+1)-S(K)*(X(K)-X(K+1)))/(S(K+1)-S(K))
06500 AX=X(K+1)+U
06600 AY=Y(K+1)+S(K+1)*U
06700 XA=X(K)
06800 1610 XB=X(K+1)
06900 IF(XA.LE.XB) GO TO 202
07000 XA=X(K+1)
07100 XB=X(K)
07200 202 YA=Y(K)
07300 YB=Y(K+1)
07400 IF(YA.LE.YB) GO TO 204
07500 YA=Y(K+1)
07600 YB=Y(K)
07700 204 CONTINUE
07800 IF(AX.GE.XA.AND.AX.LE.XB.AND.
07900 1 AY.GE.YA.AND.AY.LE.YB) GO TO 290
08000 205 K1=K1+1
08100 DO 210 K2=K+1,N
08200 K3=N+K+1-K2
08300 X(K3+1)=X(K3)
08400 Y(K3+1)=Y(K3)
08500 I1(K3+1)=I1(K3)
08600 I2(K3+1)=I2(K3)
08700 210 S(K3+1)=S(K3)
08710 IF(FLG.EQ.1..OR.FLG.EQ.2.) GO TO 280
08800 N=N+1
08900 Z0=X(K)**2-X(K+2)**2
09000 Z1=Y(K)-Y(K+2)-S(K)*(X(K)-X(K+2))-((X(K)+X(K+2))/2.
09100 1 -X(K))*(S(K)-S(K+2))
09200 Z2=X(K)**3-X(K+2)**3-3*(X(K)-X(K+2))*X(K)**2
09300 1 -1.5*(X(K)+X(K+2))*Z0+3*X(K)*Z0
09400 AZ=Z1/Z2
09500 BZ=(S(K)-S(K+2)-3*(X(K)**2-X(K+2)**2)*AZ)/
09600 1 (2*(X(K)-X(K+2)))
09700 CZ=S(K)-3*AZ*X(K)**2-2*BZ*X(K)
09800 DZ=Y(K)-AZ*X(K)**3-BZ*X(K)**2-CZ*X(K)
09900 X(K+1)=-BZ/(3.*AZ)
10000 Y(K+1)=AZ*X(K+1)**3+BZ*X(K+1)**2+CZ*X(K+1)+DZ
10100 S(K+1)=3*AZ*X(K+1)**2+2*BZ*X(K+1)+CZ
10200 IF(Y(K+1).LT.YB.AND.Y(K+1).GT.YA) GO TO 290
10300 X(K+1)=(X(K+2)+X(K))/2.
10400 Y(K+1)=(Y(K+2)+Y(K))/2.
10500 S(K+1)=0.
10510 GO TO 290
10520 280 N=N+1
10530 IF(FLG.EQ.2.) GO TO 282
10540 S(K+1)=2.*(Y(K+2)-Y(K))/(X(K+2)-X(K))
10550 GO TO 284
10560 282 S(K+1)=(Y(K+2)-Y(K))/(2.*(X(K+2)-X(K)))
10570 284 X(K+1)=(X(K+2)+X(K))/2.
10580 Y(K+1)=(Y(K+2)+Y(K))/2.
10600 290 IF(S(K).EQ.-1001.) S(K)=1001.
10700 IF(S(K+1).EQ.-1001.) S(K+1)=1001.
10800 300 CONTINUE
10900 C FIT THE POINTS WITH N-1 CURVES
11000 305 DO 400 K=1,N-1
11100 IF(I1(K).EQ.0) GO TO 310
11200 A(1,K)=0.
11300 A(2,K)=0.
11400 B(1,K)=X(K+1)-X(K)
11500 B(2,K)=Y(K+1)-Y(K)
11600 C(1,K)=X(K)
11700 C(2,K)=Y(K)
11800 GO TO 400
11850 310 IF(S(K)**2.EQ.SZ.OR.S(K+1)**2.EQ.SZ) CALL S101
11900 B(1,K)=2*(Y(K+1)-Y(K)-S(K+1)*(X(K+1)-X(K)))/
12000 1 (S(K)-S(K+1))
12100 B(2,K)=S(K)*B(1,K)
12200 A(1,K)=X(K+1)-X(K)-B(1,K)
12300 A(2,K)=Y(K+1)-Y(K)-S(K)*B(1,K)
12400 C(1,K)=X(K)
12500 C(2,K)=Y(K)
12600 400 CONTINUE
12700 C CALCULATE 512 POINTS
12800 M=N-1
12900 R=M/511.
13000 T=-R
13100 DO 500 L=1,511
13200 T=T+R
13300 KT=T+1
13400 KU=T
13500 TT=T-KU
13600 X1(L)=A(1,KT)*TT**2+B(1,KT)*TT+C(1,KT)
13700 500 Y1(L)=A(2,KT)*TT**2+B(2,KT)*TT+C(2,KT)
13800 X1(512)=X(N)
13900 Y1(512)=Y(N)
14000 600 RETURN
14100 END